Journal of Mathematical Neuroscience (2014) 4:9 O The Journal of Mathematical Neuroscience 

DOI 10.1186/2190-8567-4-9 a SpringerOpen Journal 



RESEARCH Open Access 



Identification of Criticality in Neuronal Avalanches: II. 
A Theoretical and Empirical Investigation of the Driven 

Case 

Caroline Hartley • Timothy J. Taylor • 

Istvan Z. Kiss • Simon F. Farmer • Luc Berthouze 

Received: 16 September 2013 / Accepted: 20 March 2014 / Published online: 25 April 2014 

© 2014 C. Hartley et al.; licensee Springer. This is an Open Access article distributed under the terms of 

the Creative Commons Attribution License (http://creativecommons.Org/licenses/by/2.0), which permits 

unrestricted use, distribution, and reproduction in any medium, provided the original work is properly 

cited. 

Abstract The observation of apparent power laws in neuronal systems has led to the 
suggestion that the brain is at, or close to, a critical state and may be a self-organised 
critical system. Within the framework of self-organised criticality a separation of 
timescales is thought to be crucial for the observation of power-law dynamics and 
computational models are often constructed with this property. However, this is not 
necessarily a characteristic of physiological neural networks — external input does not 
only occur when the network is at rest/a steady state. In this paper we study a simple 
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neuronal network model driven by a continuous external input (i.e. the model does 
not have an explicit separation of timescales from seeding the system only when in 
the quiescent state) and analytically tuned to operate in the region of a critical state 
(it reaches the critical regime exactly in the absence of input — the case studied in 
the companion paper to this article). The system displays avalanche dynamics in the 
form of cascades of neuronal firing separated by periods of silence. We observe par- 
tial scale-free behaviour in the distribution of avalanche size for low levels of external 
input. We analytically derive the distributions of waiting times and investigate their 
temporal behaviour in relation to different levels of external input, showing that the 
system's dynamics can exhibit partial long-range temporal correlations. We further 
show that as the system approaches the critical state by two alternative 'routes', dif- 
ferent markers of criticality (partial scale-free behaviour and long-range temporal 
correlations) are displayed. This suggests that signatures of criticality exhibited by a 
particular system in close proximity to a critical state are dependent on the region in 
parameter space at which the system (currently) resides. 

1 Introduction 

In recent years, apparent power laws (i.e. where a power law is the best model for the 
data using a model selection approach [1,2]) have been observed experimentally in 
neurophysiological data leading to the suggestion that the brain is a critical system [3, 
4]. These observations have included that of neuronal avalanches — cascades of neu- 
ronal firing recorded in vivo and in vitro whose size and duration appear to follow 
power-law distributions [5-9]. Recently it has been claimed that equivalent neuronal 
avalanche behaviour with the same power-law relationship can be identified in hu- 
man MEG (magnetoencephalography) recordings [10]. On a wider scale, fluctuations 
in oscillation amplitude in human (adult and child) EEG (electroencephalography) 
and MEG exhibit a power-law decay of the autocorrelation function of the signal — a 
property known as long-range temporal correlations (LRTCs) [4, 11-15]. These ob- 
servations and the idea that the brain is a critical system have drawn much attention as 
critical systems have been shown to exhibit optimal dynamic range and optimal infor- 
mation processing [16, 17]. Moreover, it has led to the hypothesis that brain dynam- 
ics may fit within the framework of self-organised criticality (SOC), i.e. a system that 
does not require external tuning of parameters to reach the critical state [4, 18, 19]. 

While the observation of power laws within neuronal activity may be attractive 
we must address the issue of whether (specifically) a neuronal system in the region 
of a critical state can produce this type of dynamics. Propagation of the spiking of 
neurons within a network has been interpreted within the context of percolation dy- 
namics and the theory of branching processes [20, 21]. A critical branching process 
is a process such that one active node will activate on average one other node at the 
next time step and so one can discern how this would relate to neuronal systems 
whereby the system is critical if one active neuron on average activates one other 
neuron at the next time step. A critical branching process will display power-law dy- 
namics, however, a number of assumptions underlying branching processes do not 
hold true in neurophysiological systems. Firstly, the theoretical analysis of branching 
processes relies on full-sampling of the system. Full-sampling is unlikely to occur in 
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the experimental setting and this can have a profound effect on the distribution [22] . 
Additionally, re-entrant connections invalidate the standard theory of branching pro- 
cesses [21] and so this brings into question the idea that neuronal systems can be 
modelled as critical branching processes. Moreover, the strict definition of a critical 
system is one that operates at a second-order phase transition which applies only to 
systems with infinite degrees of freedom. Therefore, we should expect a critical sys- 
tem to exhibit an exact power-law distribution in the case of infinite size but what 
should we expect if the system is finite? As neuronal systems are necessarily finite 
this is an important question in the neuroscience field but one that has yet to be fully 
addressed. Within experimental results this fact has been accounted for by the con- 
cept of finite- size effects — where a power law is observed up to a cut-off value [2, 
5, 6, 19]. This cut-off value has been suggested to coincide with the size of the sys- 
tem and distributions from networks of different sizes have been shown to exhibit 
an exact scaling relationship — a phenomenon known as finite-size scaling [2, 23]. 
However, the finite- size effect with a cut-off value at system size has been assumed 
without analytical derivation (though, see the companion paper to this article [24], as 
described below) and the questions of how a finite critical system behaves and what 
types of dynamics are possible for such a system remain open in the field. Whether 
a finite-size system should display the same signatures of criticality as the system in 
the limit of system size is not known. 

In the companion paper to this article [24] we examined a computational model of 
a finite neuronal system analytically tuned to its critical state, defined as a transcritical 
bifurcation. There we showed that the dynamics of the system, which by analogy with 
experimental neuronal avalanches could be termed avalanches (discrete cascades of 
neuronal firing), exhibited scaling which does not follow a power law but does exhibit 
partial scale-free behaviour. We were able to show that the cut-off value is approxi- 
mately the system size, as suggested experimentally by the finite-size effect, but it is 
analytically related to the lead eigenvalue of the transition matrix (the matrix of all 
possible transitions at each simulation step). This is an important observation given 
that avalanches in systems with re-entrant connections could in principle be of infi- 
nite size and yet experimental observations have suggested that neuronal avalanches 
exhibit a finite-size cut-off [2, 5]. Overall, the results suggested that finite systems at 
criticality exhibit signatures of critical systems dynamics but do not (at least in this 
instance) exhibit exact power laws as had previously been suggested. 

While the system studied in the companion paper leads us to a greater understand- 
ing of the dynamics displayed by a finite neuronal system, there is still an important 
difference between the system studied there and physiological neuronal systems. In 
the companion paper the system was seeded by setting a single neuron in the net- 
work into the active state and an avalanche was defined as the firing that occurred 
until the network returned to a stable state (the fully quiescent state). After this point 
no more firing could occur until the system was reseeded. This imposed a separa- 
tion of timescales, with all avalanches and neuronal firing occurring on a much faster 
timescale than the timescale of the 'external input' reseeding the system. Many other 
computational models have also taken this approach [18, 25, 26], with a separation of 
timescales thought to be necessary for the observation of self-organised critical dy- 
namics [23]. While a separation of timescales is likely to occur in some natural sys- 
tems such as earthquakes, where friction in the Earth's plates build up over the course 
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of years but energy is released in a matter of minutes, this is not a physiologically re- 
alistic assumption for a neuronal system. External input (be it from the environment 
or other areas of the nervous system) will not arrive only once the neuronal popu- 
lation has returned to a set state. Before physiological recordings can be interpreted 
within the field of critical systems we must address the question of the types of dy- 
namics that should be expected by not only a finite-size system but also a system that 
is driven by a physiologically realistic external input. Can a finite-size system without 
an explicit separation of timescales in the region of a critical regime exhibit markers 
of criticality? How might the external input to the system affect these markers? 

Previous authors examining computational neuronal networks with continuous 
driving (i.e. no explicit separation of timescales) have observed power-law dynam- 
ics [16, 27-29]. In particular, Kinouchi and Copelli [16] and Larremore et al. [29] 
analytically determined the parameters required such that the model they studied was 
at criticality and displayed peak dynamic range, in fully connected networks and 
networks with a range of topologies, respectively. However, these authors did not ex- 
plicitly examine the firing dynamics of the system in the region of the critical regime, 
concentrating on average activity levels. In a SOC system such as the sand-pile model 
[18] the waiting times (periods of inactivity between avalanches) have been shown to 
follow an exponential distribution [30]. However, these waiting times are related to 
the reseeding of the system — sand is added to cells chosen at random and the next 
avalanche begins when a cell exceeds the threshold. In contrast, recent experimental 
work has shown that waiting times between neuronal avalanches in cultures have a 
distribution with two trends — a (short) initial power-law region thought to relate to 
neuronal up-states and a bump in the distribution at longer waiting times thought to 
relate to neuronal down states [31]. Could this difference in these waiting time dis- 
tributions (between the SOC sand-pile model and the neuronal avalanches in culture) 
be explained by the fact that physiological neuronal systems do not have a separation 
of timescales? 

As described above, another signature of criticality that has been reported in neural 
systems is the presence of LRTCs. In the majority of cases they have been observed in 
large scale neuronal signals such as human brain oscillations. Recent endeavours have 
been made to link these observations of scale-free behaviour on large scales with neu- 
ronal avalanches [32, 33]. Poil et al. demonstrated in a computational neuronal net- 
work that power-law distributed avalanches and LRTCs in oscillations emerge con- 
currently. In addition, LRTCs have also been detected in the waiting times of bursts 
of activity in cultures [34] and the discontinuous burst activity recorded in the EEG 
of extremely preterm human neonates [35]. Thus, LRTCs have been demonstrated 
in discrete neuronal activity yet they have not been examined in the waiting times 
of neuronal avalanches themselves. While LRTCs in avalanche activity would not be 
possible in a seeded computational system (where the activity is initiated 'by hand' 
and there is no memory within the system's dynamics) it is conceivable that a driven 
system, which is more akin to physiological networks which can display LRTCs, 
might display this type of dynamics in the waiting times of neuronal avalanches. 

In summary, in this paper we aim to address the following questions: 

1 . Assuming that the brain, or population of neurons under study, operates in the re- 
gion of a critical regime can it be expected to display power-law statistics given 
that it is a finite-size system? If not what distribution should we expect? As dis- 
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cussed, this question was also addressed in the companion paper [24], where we 
studied a system without an external input. However, here we specifically consider 
this question in the context of a driven system. 

2. Can we expect a finite- size neuronal system in the region of a critical regime to 
exhibit other markers of criticality, and specifically the presence of LRTCs? Does 
the presence of LRTCs relate to that of power-law distributions? As described 
above, LRTCs have been observed in neurophysiological data sets. However, a full 
theoretical examination of how LRTCs may relate to other markers of criticality 
in neuronal systems is lacking. 

3. How are signatures of criticality (power-law distributions and LRTCs) affected by 
proximity to the critical regime? One might assume that a system which is closer 
to a critical regime may exhibit signatures of criticality, whereas a system that is 
further from the critical regime will not. Importantly, our analysis shows that this 
assumption is in fact not (always) true. 

Although these questions are particularly applicable and novel to the field of neu- 
roscience, it should be noted that similar questions have been the subject of much 
research in the field of statistical physics; see [36] for one review. Moreover, Marko- 
vian neural models with saturating firing functions have been suggested previously 
to fall within the universality class of directed percolation when analysed at the con- 
tinuum and thermodynamic limit [37]. However, it is important to make the distinc- 
tion between criticality in the statistical mechanics sense (i.e. defined in terms of 
a second-order phase transition in a system with infinite degrees of freedom) and 
criticality in the mathematical sense (i.e. defined in terms of a bifurcation in a low- 
dimensional mean field model) such as used in our work. This paper will show some 
shared phenomenology although it should be clear that properties, in particular, uni- 
versal properties, associated with some classes of critical phenomena in the statistical 
mechanics sense should not necessarily be expected from either our or other related 
neuroscience models. 

In this paper, as in the companion paper, we examine a purely excitatory (in terms 
of synaptic transmission — see Discussion) stochastic neuronal model. As in the com- 
panion paper, a number of assumptions are made to simplify the model with the 
outcome that it is analytically tractable and therefore can be tuned to operate in the 
region of a critical regime. This approach is taken as it allows direct exploration of 
the above questions, which would not be possible with a more complex system. We 
begin by examining the distributions of avalanche size and duration, investigating 
the presence of scale-free behaviour. We also show that as the system approaches 
the theoretical critical regime by decreasing the external input, there is a change in 
the distributions of avalanche characteristics with the appearance of partial scale-free 
behaviour in avalanche size. It is important to note that the definition of avalanches 
strongly depends on the choice of binning method. In the literature different defini- 
tions of avalanches are used in models with seeded systems and with systems where 
the dynamics is continuous (including physiological recordings). We will return to 
this in the Discussion. 

Unlike in the companion paper where the system was seeded after each avalanche, 
the system studied here was driven by a continuous external input. This allowed us 
to additionally assess the waiting times, which are intrinsic to the system, and we 
were able to analytically derive the distribution of waiting times. We then investi- 
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gated the presence of partial LRTCs in the empirically derived waiting times. Finally, 
we showed that as the system size increases (and the system approaches the theo- 
retical critical regime from a different route) the range over which the correlations 
extend also increases. Overall we find that the system displays different signatures of 
criticality depending on the region of the parameter space around the critical regime. 



2 The Model 

In this paper, as in the companion paper, we study a stochastic model based on that 
of Benayoun et al. [38], an extended version of the previously introduced stochastic 
rate model [37, 39]. Though greatly simplified from a physiological neural network, 
the model is chosen as it is analytically tractable and thus enables direct derivation 
of the parameters such that there is a critical regime. With this approach it is there- 
fore possible to assess the dynamics of a neuronal system in the region of (or at) 
a critical regime. While Benayoun et al. considered a network with both excitatory 
and inhibitory connections, we simplify the system further (as in the companion pa- 
per), considering a network with purely excitatory synaptic connections. As will be 
discussed later, this type of network can be set within the context of early brain de- 
velopment. 

We consider a system of TV fully connected neurons, with each neuron in one of 
two states — active (A) or quiescent (Q). For a small time step dt —> 0 the probability 
of transition for a neuron between the two states is given by 

P(Q -> A, in time dt) = f(si(t))dt, 
P(A -> Q, in time dt) = adt, 

where Si (t) = • ^jfcij it) + hi (t) is the input to neuron i , f is an activation function, 
hi it) is the external input to neuron i, Wij is the connection strength from neuron i 
to neuron j and a jit) = 1 if neuron j is active at time t and zero otherwise. Finally, 
a is a constant rate at which neurons change from the active to inactive (quiescent) 
state. 

For analytical tractability and characterisation of the critical state we make the 
following additional simplifications: 

1. The synaptic connection strengths are the same for all connections with Wij = 
w > 0. 

2. The external input is constant to all neurons and at all simulation steps so that 
hiit) = h>0. 

3. The activation function is linear with f(x)=x. 

While the first and third assumptions are the same as in the companion paper, we 
make the additional assumption of constant positive external input here as opposed to 
the companion paper where we examined the system with no external input (h = 0). 
As the network is fully connected, and the system is closed so that A + Q = N 
(where A is the number of active neurons and Q is the number of quiescent neurons), 
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the system can be described by the mean field equation: 

dA (wA \ 

- = (- +h y N - A) -aA. 

As stated in the companion paper, we can use this equation to analyse the stability of 
the system about the fixed point and determine the parameters for which the system 
is at the threshold of stability, i.e. when the fixed point is critical. This threshold 
occurs when the eigenvalue (A) of the fixed point is zero, which can alternatively 
be stated, borrowing terms from the epidemiology literature, as Ro = 1 (the basic 
reproductive ratio). Moreover, this is also equivalent to a branching parameter of one. 
In the companion paper it was shown that with h = 0, Ro = ^ and so for Rq = ^ = 
1 =>> a = w the system is critical. 

Here we study the system in the presence of a positive external input, h > 0. In 
this case the fixed point of the system is given by 

— — A 2 + wA + h(N - A) - a A = 0, 
N 

and the eigenvalue at the fixed point is 

w 

A = —2 — A + w — h — a. 
N 

For a fixed point to be critical we require that both these equations be satisfied. 
However, solving them simultaneously we find that there are no real roots when 
w, h, N > 0. This implies that there is no parameter region such that the system (with 
this activation function and positive external input) has a critical fixed point. How- 
ever, considering again the case with no external input (h = 0) for which the critical 
state occurred with parameters a = w, if this system is driven by a 'sufficiently low' 
level of external input it should still be within the region of the critical state. There 
has been some suggestion that the brain is not directly at a critical point but is in fact 
just very close to the critical regime and it has been speculated that the brain may 
actually be slightly supercritical [32]. Additionally, it been shown that a computa- 
tional model of neuronal avalanches which follows a SOC approach [25] is actually 
a system that 'hovers' close to the critical state [23]. Therefore, the question of how 
a finite driven system behaves in the region of a critical regime is pertinent to the 
neuroscience field. 

An additional motivation for considering a non-zero external input is the dynamic 
range (A) of the system. Larremore et al. [29] describe the dynamic range as "the 
range of stimuli over which there is significant variation in the collective response 
of the network". Kinouchi and Copelli [16] examine dynamic range in models of 
networks with uniform connectivity operating with discrete time dynamics where 
multiple firings can occur within each time step. They found that the dynamic range 
was maximised when the local branching ratio was equal to one. Larremore et al. 
[29] consider a version of this model but with the introduction of heterogeneity in 
connections, showing that it is the lead eigenvalue, A, of the connectivity matrix that 
governs the dynamic range and that the dynamic range is maximised when A = 1 . 
In Appendix A we provide an analytic calculation for the dynamic range of our con- 
tinuous model. Analogous to the results described above [16, 29] the dynamic range 
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Rq 

Fig. 1 Rq versus A. Plot of Rq versus A (the dynamic range — see Appendix A) for both the analytic 
result (black line) and simulations (red dots). For the comparable simulated result we average 10,000 
realisations that have run until time t = 200. To obtain a reasonable spread of h we used the conjugation 
of the intervals [0 : 0.002 : 0.2] and [0.4 : 0.2 : 18] 

is maximised when Ro = 1 (w/a = 1 when h = 0). This is illustrated in Fig. 1 where 
results from simulations are compared with the analytic solution. It is important to 
emphasise that the parameterisation of the dynamic range is in terms of the value Ro 
calculated for networks when there is no external input. When this parameterisation 
is such that Ro = 1 =>> a = w (and therefore when the system is tuned to the critical 
state) external input to the system will give rise to dynamics for which the dynamic 
range is maximised. This point will be considered further in the Discussion. 

Throughout this study we will examine the system in the presence of an external 
input of h = l/N or less. This level of the external input is equivalent to setting 
a single neuron to the active state and so corresponds to seeding the system in the 
zero input case. We therefore deem this level of the external input to be sufficiently 
low such that we would expect the system to remain within the region of the critical 
regime. As in the companion paper we set w = a = 1 . With these parameters and 
with positive external input we find that the fixed point of the system is given by 

hN lN 2 h 2 ~7 
A = ± J + N 2 h 

2 V 4 

(as the solution A < 0 is not physical, it is ignored, leaving a single positive fixed 
point for positive h), and the eigenvalue of this fixed point is given by 

X = -y/h 2 + 4h. 

With lower levels of external input the system approaches the critical regime (see 
Fig. 2). Note that this approach is in fact from a slightly subcritical state given these 
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Fig. 2 The eigenvalue of the system compared with the level of external input and system size, a With 
w = a the eigenvalue decreases with lower levels of external input h, with the system reaching the critical 
regime in the absence of external input (A, = 0 i.e. Rq = 1 , the case studied in the companion paper), b With 
h = 1 /N and w = a the eigenvalue of the fixed point X -> 0 as TV -> 00. Thus, the system approaches the 
critical state as the system size increases 

values of a and w and with positive external input. Under these conditions it is not 
possible to consider an approach from a supercritical regime with a positive eigen- 
value. 

As described above, we (initially) set h = 1 /N. With this level of external input 



As N -> 00, X -> 0 (see Fig. 2). Thus, for this level of the external input (h = l/N), 
as the system size (N) increases the system approaches the critical state (as the system 
reaches the critical state exactly when the eigenvalue A = 0). We will examine the 
effect on the dynamics of decreasing the external input, thereby allowing the system 
to approach the critical regime. We will also investigate an alternative route to the 
critical regime by increasing the system size in systems with a constant (overall) 
level of external input. 

2.1 Model Simulations and Burst Analysis 

As in the companion paper and in Benayoun et al. [38], simulations of the network 
dynamics were carried out using the Gillespie algorithm for stochastic simulations 
[40] . Briefly, at each step in the simulation 

• The total transition rate r for all the neurons within the network is calculated, with 
r = r aq + r qa where r aq is the total rate of active -> quiescent transitions and is 
given by r aq =aA and r qa is the total rate of all quiescent -> active transitions 
which is given by r qa = f(si)(N - A). 

• The time to the next transition dt is selected at random from an exponential distri- 
bution of rate r. 

• The type of transition is selected by generating a random number n e [0,1]. If 




and the eigenvalue of the fixed point is given by 
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200 400 600 800 1000 
Time (ms) 



Fig. 3 Raster plots of the network dynamics for different levels of the external input. Neuronal firing 
across 1 second of a simulation with an external input of h = 1 /N (a), h = 0. 1 /N (b) and h = 0.01 /N (c). 
For all three simulations TV = 800, a = w = 1 and the red line indicates the rate of firing in 1 ms bins. As 
can be observed, as the level of the external input decreases the firing rate decreases and the time between 
avalanches increases 



n < ^ then a randomly chosen active neuron becomes quiescent, otherwise a 
(randomly chosen) quiescent neuron switches to the active state. 

At each step in the simulation a single neuron makes a transition, though the rate at 
which transitions occur changes and so the simulation step changes. If the network 
is in a fully quiescent state (Q = N) then, with positive external input, r aq = 0 but 
r qa = hN and consequently there will necessarily be a transition of a neuron from 
the quiescent to the active state. Similarly, when the network is in the fully active 
state (A = N) r qa = 0 but r aq = otN and so there will necessarily be a transition 
of a randomly chosen neuron from the active to the quiescent state. From all other 
starting points transitions from the active to the quiescent or from the quiescent to the 
active state are possible. Thus, from all network states one neuron will change state. 
This is unlike the companion paper where with no external input the network must 
be seeded when in the fully quiescent state. Instead in this case network dynamics is 
continuous (i.e. no re-seeding is required) and are of finite length only in-so-far as 
they are restricted by simulation lengths. 

We defined a neuron as firing at the first time step at which the neuron switches 
from the quiescent to the active state. Figure 3 shows raster plots of network firings 
for the first 1 second of simulations with three different levels of the external input. 
As was described above, unlike in the companion paper where there was no external 
input, the dynamics continues even if the system reaches the fully quiescent state. In- 
terestingly, we can also notice that the network dynamics appears to exhibit burst like 
behaviour, with periods with high neuronal firing interspersed with periods without 
network firing. It is important to realise that these bursts are intrinsic to the system 
and are not directly related to the dynamics of the external input (the input is constant 
to all neurons in each of the simulations) nor due to a saturation of the network — the 
bursts themselves consist of different numbers of neuronal firing. In all three cases 
the parameters are set to the critical state (with no external input). With lower levels 
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of the external input the system approaches the critical regime and we see that the 
bursts become further apart and more distinct. We will examine the distributions of 
this dynamics below. (See also Appendix B where we examine driving the system 
from subcritical and supercritical states.) 

This burst dynamics is analogous to the neuronal avalanches observed experimen- 
tally in that they are discrete cascades of firing. Neuronal avalanches observed ex- 
perimentally in physiological networks are so called because they have sizes which 
are distributed according to a power-law and while the size distribution of the burst 
activity in this network has yet to be presented we will refer to the activity throughout 
the rest of this paper as avalanches due to their discrete burst behaviour. To determine 
the distribution of the avalanches we divided the activity into individual avalanches 
using the approach of Benayoun et al. [38]. This method divides consecutive neu- 
ronal spiking between any two neurons within the network into separate avalanches 
if the time difference between the spikes is greater than the average difference (8t) 
between consecutive spikes within the simulation. This approach (referred to later 
in the text as the binning method) is similar to the method used to define neuronal 
avalanches within physiological data [5, 6] — though the choice of binning method 
will be discussed later in the paper. It is important to note that this binning approach 
used to define avalanches was not used in the companion paper, where an avalanche 
was defined as all firing that occurred before the network reached the fully quiescent 
state and was reseeded. This has been used as a standard classification for discon- 
tinuous data, stemming from the sand-pile model of criticality [18]. However, as the 
firing dynamics here continues for the entire simulation it was instead appropriate to 
use an approach that had been used previously for continuous dynamics. 

Throughout the remainder of this paper we examine characteristics of these 
avalanches: namely the size and duration of avalanches as well as the inter- avalanche 
intervals (IAIs). The size of an avalanche is defined (in the standard way) as the num- 
ber of firings within the avalanche. If a single neuron fires more than once within a 
single avalanche it is also counted more than once. The duration of an avalanche is 
defined as the time between the start of the avalanche (the first neuron firing) and the 
end of the avalanche. Note that if the avalanche consists of a single neuron firing then 
the duration of the avalanche is 0 (and the size of the avalanche is 1). Similarly, an 
IAI is defined as the time between the end of one avalanche and the start of the next 
avalanche, i.e. the waiting time between avalanches. Note that the minimum IAI is 
bounded below by 8t as a separation between two consecutive spikes of greater than 
8t defines separate avalanches. 

2.2 Distributions of Avalanche Size and Duration 

Figure 4 shows the distributions of avalanche size and duration from example sim- 
ulations for the three different levels of external input investigated. (Compare also 
with Fig. 3 of the companion paper [24] which shows the avalanche size distribution 
in the case without external input.) With lower levels of external input the system 
approaches the critical regime and the distributions of avalanche size appear scale- 
free across a range of scales, with linearity on double logarithmic axis — the linear fits 
are indicated on the plots. The distribution approaches the distribution found in the 
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Fig. 4 Distributions of avalanche size and duration for varying levels of the external input. The distribu- 
tions of avalanche size (a, b, c) and duration (d, e, f) from simulations with h = 1/N (a, d), h = 0.1 /N 
(b, e) and h = 0.01/ TV (c, f). As the level of the external input is decreased the system approaches the 
critical regime. For all simulations TV = 800, a = w = 1 and the distributions are pooled from 10 simula- 
tions each of length 10 4 seconds. The red lines indicate linear fits on the double logarithmic scale (where 
appropriate), i.e. fitted power laws with exponents of 1.68 (b), 1.48 (c), 1.88 (e) and 1.64 (f) 



companion paper for the system exactly tuned to the critical state. With h = 0.01 /N 
(i.e. the lowest level of eternal input) the exponent of the fitted power law is approxi- 
mately 1.5, see Fig. 4(c), which is consistent with experimentally observed neuronal 
avalanche sizes [5, 6]. However, for higher levels of the external input this scale- 
freeness of the distribution is lost which coincides with moving away from the crit- 
ical regime. In the case of avalanche duration a similar relationship with the critical 
regime is seen with a scale-free portion in the middle ranges of the distribution (be- 
tween approximately 2 and 50 ms) with lower levels of external input. Thus, for 
lower levels of external input, when the system approaches the critical regime, the 
distributions, in particular the distribution of avalanche size, exhibit partial scale-free 
behaviour. 

It is worth considering what leads to the changes seen in the distributions as the 
level of external input is varied. As stated, as the level of external input decreases, the 
system approaches the critical regime and so it is perhaps not surprising that signa- 
tures of criticality (i.e. scale-free behaviour) emerge in the distribution of avalanche 
size as the external input is lowered. Examining the raster plots of firing for the dif- 
ferent levels of external input, see Fig. 3, we see that for lower levels the avalanches 
are further apart and more distinct. While the external input itself is continuous, 
at the lower levels of external input there is a separation of timescales, where one 
avalanche always finishes well before the next avalanche begins. The distribution 
therefore appears to follow similar characteristics to a system with a built in separa- 
tion of timescales and we confirm that the distribution is similar to that found in the 
companion paper (in which the model had an explicit separation of timescales, i.e. the 
system was only seeded once it had reached the quiescent state) where an exponent 
close to 1.5 was also observed for the distribution of avalanche size. As the external 
input is increased there are no longer such distinct periods between avalanches. This 
leads to a superposition effect, with the next (actual) avalanche starting before the 
previous avalanche has finished (i.e. a new network cascade is initiated before the 
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previous one has finished). This leads to these 'avalanches' being defined using the 
binning approach as a single avalanche (see Discussion). The scale-free behaviour in 
the distributions of avalanche size and duration is therefore lost. 

2.3 Theoretical Derivation of the Distribution of the IAIs and Comparison with 
Simulated Data 

The temporal patterning of activity within networks of neurons has long been investi- 
gated as a property of key importance, with neural rate and temporal coding suggested 
as potential substrates for information propagation. While it remains to be fully deter- 
mined how different neuronal firing properties may lead to information transfer this 
suggests that in addition to the distribution of avalanches sizes the intervals between 
avalanches need to be considered as a functional entity in their own right. As well as 
determining the IAI distributions through simulations we found it is possible to de- 
rive the theoretical distribution. In this section we derive this theoretical distribution 
and compare it with results from simulations. 

We begin by noting that a single IAI is a period during which there is no neuronal 
firing, i.e. neurons can only be switching from the active to the quiescent states or 
an IAI may be a period with a single quiescent to active transition which is preceded 
by another quiescent to active transition. We wish to derive the distribution of these 
periods. Let us initially ignore the fact that there is a minimum duration (8t) of an 
IAI and first consider the distribution of all consecutive active to quiescent transitions 
(we will return to the distribution of single quiescent to active transitions later). 

2.3.1 Distribution of Consecutive Active to Quiescent Transitions 

Let No be the number of active neurons at a time point in the simulation. After a 
single simulation step the number of active neurons will be No + 1 or No — 1 , as at 
every simulation step only one neuron makes a transition. Let qt be the probability 
that an active neuron goes back to the quiescent state given that there are i active 
neurons. Note that from the transition rates: 

ai aiN 
qi ~ ((w/N)i+h)(N -i)+ai ~ (wi + hN)(N — i) + aiN ' 

Starting with No active neurons the probability that there are No — I active neurons 
after a single simulation step is q^ 0 and the probability that there are No + 1 active 
neurons is 1 — q^ 0 . Given these probabilities we can construct a probability tree, 
shown in Fig. 5, particularly concentrating on the portion of the tree corresponding to 
active to quiescent transitions, i.e. those transitions that form a period of consecutive 
active to quiescent transitions. (Note that this probability tree focuses on different 
aspects of the model to that of the probability tree in the companion paper.) 

From this tree approach we can calculate that the probability of exactly k consecu- 
tive active to quiescent transitions (note that to consist of exactly k active to quiescent 
transitions the transition sequence must be ended by a quiescent to active transition): 

k-l 

P(IAI k ) = p(N 0 , k) = (1 - q No - k ) Y\ q No -j . (1) 

7=0 
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Fig. 5 Probability tree of 
consecutive active to quiescent 
transitions. Starting from a state 
with Nq active neurons the 
probability tree diagram 
indicates the possible transitions 
from each state specifically 
concentrating on the active to 
quiescent transitions. The 
probability qi of each transition 
is as indicated in the main text 
and is dependent on the number 
of active neurons, i 




The duration of this period of consecutive active to quiescent transitions is given 
by the sum of the times for each of these k transitions (plus the time for the quies- 
cent to active transition). As the Gillespie algorithm is used for simulations, at each 
simulation step the time to the next transition is drawn randomly from an exponential 
distribution with rate r (see above), where r is dependent on the number of active 
neurons and so changes at each simulation step. The duration of consecutive active to 
quiescent transitions is therefore the sum of exponentially distributed variables drawn 
from distributions of different rates, i.e. the distribution of consecutive active to qui- 
escent transitions is a hypoexponential. Thus, the duration distribution, f(x,No,k), 
of consecutive active to quiescent transitions of length x, consisting of k transitions, 
ending with an additional quiescent to active transition and starting from No active 
neurons is [41]: 

/(x^o,t) = i;r^-^( ft - ^ \ (2) 

7=0 \i=o,i^j rN °- 1 ~ rN o-J/ 

when r# 0 -i ^ r N 0 -j and where r m is the total transition rate for all neurons within 
a network with m active neurons and is the rate of the exponential distribution from 
which the time to the next transition is randomly drawn. This equation holds provided 
that r# 0 _j 7^ rN Q -j Vz, j. If this is not the case and there exists A, B e {I, .. . , N] 
such that (a + w — h)/w = A -\- B ^ ta = rs then we instead use the more general 
form — assuming there are a distinct rates, which we label P\, Pi, - - , Pa mat occur 

c\ , C2, . . . , c a times, respectively (i.e. c\ + C2 H h c a = k + 1), then the duration 

distribution is given by [42] : 

f{x,NQ,k) = By > — , (3) 

tlh (<* -/)!(/ -1)! 
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Fig. 6 Probability tree of all possible transitions in a network of size N = 3. Simulations start from a 
state with no active neurons: Nq = 0. The diagram shows the possible transitions at each step, along with 
the probability of making that transition. The probabilities are as defined in the main text with qi being 
the probability of a neuron switching from the active to the quiescent state given i initially active neurons. 
Dotted lines indicate transitions that are already shown elsewhere in the tree and so the tree shown here 
completely describes all possible transitions in a network of this size 

where 

B = Y[P C j j and 0^(0 = — j (Pj+tr Cj . 

Whilst this involves higher-order derivatives a closed-form solution is provided by 
Amari and Misra [43]. 

From Eq. 1 we know the probability of k consecutive active to quiescent tran- 
sitions. This equation holds true for any k up to k = No, which is the maximum 
number of consecutive active to quiescent transitions as the fully quiescent state is 
then reached. Therefore, the distribution, F(x, No), of consecutive active to quies- 
cent transitions of duration x starting with No active neurons but consisting of any 
number of transitions is a weighted sum of hypoexponentials: 

N 0 

F(x, N 0 ) = /(*> N 0 , k)p(N 0 , k). (4) 
k=\ 

2.3.2 Probability Distribution of the Initial Number of Active Neurons 

Finally, to calculate the full probability distribution of consecutive active to quiescent 
transitions for a network of set system size, N, we must combine Eq. 4 with the prob- 
ability of the initial number of active neurons being equal to No e {1,2, . . . , N] (note 
that TVo = 0 is not considered as the next transition will necessarily be an activation). 
To determine this probability, first let us consider the simple case of N = 3. We as- 
sume that the simulation starts from a state with no active neurons. Figure 6 shows all 
possible transitions between the number of active neurons in a network of this size. 
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From this probability tree the probabilities, P(i), of the number of active neurons 
being equal to i, where i e {0, 1, 2, 3} are given by 

P(0)= qi P(l), 

P(l) = P(0) + q 2 P(2), 

(5) 

P(2) = (l-$i)P(l) + P(3), 
P(3) = (1-^ 2 )P(2) 

(assuming a steady state has been reached such that the probabilities are time inde- 
pendent). Rearranging and substituting to write the equations in terms of P(l): 

P(0) = qi P(l), 
(1 -q\) 

P(2) = - — P(l), (6) 

qi 

(\-qi)(\-q 2 ) 
P(3) = — — P(l). 

42 

Furthermore, the sum of all the probabilities must equal 1 and so 

(1 -qi) (1 -qi)(l -qi)\ 
qi + l + — + — — P(l) = 1. (7) 

qi qi ) 

Therefore, 

P(l) = — . (8) 

<7i<72 + qi + (1 - ^l) + (1 - 4i)(l - ^2) 

By substituting this value back into the set of equations (Eq. 6) the probabilities for 
the full system can be calculated. 

2.3.3 Generalisation to a System of Any Size N 

From considering this simple example we can extend this to derive the probabilities 
of the number of active neurons for a system of any size TV. Firstly, as in Eq. 4 the 
probabilities can be written as (again assuming a steady state): 

P(0) = 4iP(D, 

P(1) = P(0) + 4 2 P(2), 

P(2) = (l- qi )P(l) + q 3 P(3), 



P(k) = (1 - q k -l)P(k - 1) + qk+lPik + 1), 

P(N - 1) = (1 - q N -2)P(N - 2) + q N P(N), 
P(N) = (l-q N - l )P(N-l), 



(9) 
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where = 1 but it will remain in the equations so as to aid notation. Rearranging 
gives 



(1 -q\) 
P(2) = - — P(l), 

qi 



and by induction: 



P(k + 1) = — (P(k) - (1 - q k -i)P(k - 1)) 

qk+\ 

a-qi)(l-q 2 )---(l-qk) 
= - — — — P(l). (10) 

qiqz • • • qk+\ 

Summing all the probabilities and setting this equal to 1 : 

= q 2 q3 ' ' ' qN/(q\qi ■ ■ ■ qN + <72<?3 --qN 

+ (1 - qi)q 3 ■ • ■ q N + • ■ • + (1 - qi)(l - qi) ■ • • (1 - q N )). (11) 

Having determined these probabilities we then need to take into account the fact 
that the consecutive active to quiescent transitions must be preceded by a quiescent to 
active transition, i.e. they must be preceded by a neuron firing (otherwise they would 
be a chain of k + 1 consecutive active to quiescent transitions and so included else- 
where in the distribution). We are therefore only interested in the probability Pa(No) 
of the number of active neurons being equal to No given that a quiescent to active 
transition has just occurred. Considering again the probability tree, Fig. 6, we find 
that these probabilities, are given by 

Pa(0)=0, 
P A (1) = P(0), 



(12) 

P A (k) = (l-q k - 1 )P(k-l), 



P A (N) = (l-q N -i)P(N-l), 

where we make use of the previously defined probabilities P (k) . From these probabil- 
ities the full probability distribution of the duration of consecutive active to quiescent 
transitions can be calculated. As was shown above, for a set initial number of active 
neurons TVo, the probability distribution of consecutive active to quiescent transitions 
is given by a weighted sum of hypoexponentials; see Eq. 4. This can then be further 
weighted by the probability Pa(No) that the initial (at the start of the sequence of 
transitions) number of active neurons is equal to TVo and the previous transition was 
quiescent to active. Thus, the overall probability distribution of consecutive active to 
quiescent transitions is given by 

i=0 \ m =\ I 
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Fig. 7 Theoretical and simulated distributions of periods of consecutive active to quiescent transitions. 
The simulated distribution (black) is compared with the theoretically derived probability density function 
(shown in red; see Eq. 13). For both distributions a = l,w = l,N = 50 and h = 1/N. The mean difference 
between consecutive spikes (St) within the simulation (green) is used to define avalanches through the 
binning approach described in the main text. Thus, the portion of the distribution for which the length 
of active to quiescent transitions are greater than this average time between consecutive spikes form the 
distribution of IAIs when combined with the distribution of single quiescent to active transitions 

To confirm that this theoretically derived distribution compares with results from 
simulations, we determined the distribution of the lengths of periods of any consec- 
utive active to quiescent transitions from simulations. Figure 7 shows the distribu- 
tion of consecutive active to quiescent transitions from a simulation with N = 50 
compared with the theoretical distribution; we can see that there is good agreement 
between the two. 

2.3.4 Distribution of Single Quiescent to Active Transitions 

As was described above, a single period in between two neurons firing can also be an 
IAI (providing that the duration is longer than the average time between spikes as ac- 
counted for below, note that only single periods are considered as consecutive periods 
necessarily include neurons switching to the active state and so cannot form part of 
an IAI). Thus, the distribution of IAIs should also take into account the distribution 
of single quiescent to active transitions. As we make use of the Gillespie algorithm, 
the duration distribution of these single transitions is an exponential with rate given 
by the total transition rate, which is dependent on the number of active neurons, No. 
This is then weighted by the probability of a quiescent to active transition given No 
active neurons (i.e. by 1 — q^ Q ) and additionally weighted by the probability of start- 
ing with No active neurons following a quiescent to active transition as calculated 
above. Thus, the probability distribution of single quiescent to active transitions of 
length x is given by 



i=0 

Figure 8(a) shows the simulated distribution of single quiescent to active transitions 
compared with the theoretical distribution. 

2.3.5 The IAI Distribution 

As discussed above, the IAI distribution combines these two distributions — the dis- 
tribution of consecutive active to quiescent transitions and the distribution of single 




(14) 
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Fig. 8 Theoretical distributions of single quiescent to active transitions and the combined distributions. 
Theoretical (red) and simulated (black) probability distributions for a single quiescent to active transitions 
and b this distribution combined with the distribution of consecutive active to quiescent transitions (see 
Fig. 7). In both cases a = 1, w = 1, N = 50 and h = 1/N. The green line indicates the average time 
between consecutive spikes (St) within the simulations. Thresholding the combined distribution (b) at this 
level determines the IAI distribution 

quiescent to active transitions. This combined distribution, along with simulated val- 
ues, is shown in Fig. 8(b). As was described above, avalanches are defined from the 
network firing pattern as consecutive spikes where the time difference between them 
is no greater than the average time difference between consecutive spikes, 8t, within 
the network. Thus, the minimum IAI is bounded below by 8t and all consecutive 
active to quiescent transitions or single quiescent to active transitions whose total 
duration is greater than 8t will be an IAI. Thresholding the combined distribution 
at 8t determines the IAI distribution. Figure 9(a) shows theoretical and simulated 
IAI distributions displayed on a double logarithmic scale. Despite the fact that the 
distribution is not a power law (theoretically we know that it is a weighted sum of 
hypoexponentials), it appears scale-free over a range of scales on this double loga- 
rithmic scaling. As we will see below, the distribution can also pass statistical tests for 
power-law distributions, indicative of the partial scale-free behaviour of the system 
close to the critical regime. 

Figure 9 also shows the theoretical and simulated distributions for lower levels 
of external input. With lower levels of external input (as the system approaches the 
critical regime) the average IAI increases and the distribution changes, no longer 
exhibiting scale-free behaviour. Even if we consider the same scale for all levels of 
the external input (IAIs in the region of 0.05-5 ms) then only with h = 1/N is the 
distribution scale-free. Indeed, for the lowest level of h = 0.01 /N the distribution is 
in fact well fit by an exponential, in this case y = 0.028^ _0 01x as seen in Fig. 10, 
indicating the loss of the scale-free behaviour in the distribution. Thus, scale-free 
behaviour in the case of the IAI distribution does not increase with proximity to the 
critical regime. 

As an aside, note that due to the product in the hypoexponential (see Eq. 2) de- 
termination of the probabilities for large N can become computationally intractable. 
For simulations larger than with N = 50 we therefore only determined the theoretical 
distribution up to a set level of the number of active neurons. We set the threshold 




Springer 



Page 20 of 42 



C. Hartley et al. 



(a) (b) (c) 




d.05' 1 o.5o' 1 5.bo' *~ b:5' 1 5:0 ' ' 5(3.0' 6 '2'0 '100' )doc 

IAI IAI IAI 



Fig. 9 Distribution of IAIs for varying levels of the external input. The theoretical (red) and simulated 
(black) IAI distributions with h = l/N (a), h = 0.\/N (b) and h = 0.01/ TV (c). These distributions were 
for TV = 800 with a = w = 1, with a simulation length of 10 4 seconds. The distributions from the simulated 
data are pooled from 10 simulations. The theoretical distributions were calculated up to the level of active 
neurons which occur with a cumulative probability of 0.9 (see main text). The blue line in a indicates a 
linear fit, i.e. a fitted power law with an exponent of 2.71 
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Fig. 10 The IAI distribution with h = 0.01/ N is well fit by an exponential distribution. The theoretical 
IAI distribution at the lowest level of the external input (shown in red; see also Fig. 9(c)) compared with 
the fitted exponential distribution (black dashed). The exponential is given by y = 0.028e-° Mx . This 
indicates that as the external input is decreased and the system approaches the critical regime the IAI 
distribution loses the scale-free behaviour seen at higher levels of the external input and is dominated by 
an exponential 

level of the number of active neurons according to the probability distribution of 
starting from a particular number of active neurons (calculating the cumulative prob- 
ability from zero active neurons), and sufficiently low so that the calculations were 
computationally viable. However, the theoretical distributions calculated using this 
threshold are still a good fit to the simulated data — see Fig. 9. 

2.3.6 Distributions of Avalanche Size and Duration 

As we have shown, the theoretical distribution of IAIs can be calculated by assessing 
periods of consecutive active to quiescent transitions and single quiescent to active 
transitions. It is also possible to derive the distribution of consecutive quiescent to 
active transitions. However, if a period of active to quiescent transitions (a period 
without firing) has a duration less than the average time difference between two spikes 
then this interval does not separate an avalanche into two. Therefore, the distributions 
of number and length of consecutive quiescent to active transitions does not describe 
the distributions of avalanche size and duration — these distributions can also contain 
periods of active to quiescent transitions within two or more periods of quiescent to 
active transitions. Note that a period of active to quiescent transitions having a length 
less than the average difference between consecutive spikes is not dependent on the 
number of active to quiescent transitions within the interval, as the length of each 
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transition is drawn at random from an exponential distribution. It was therefore not 
possible for us to determine a theoretical distribution of avalanche size and duration 
using this approach. 

2.4 Statistical Comparison with a Power-Law Distribution 

The influential paper by Clauset et al. [1] developed a model selection based method- 
ology to determine whether empirical data is likely to be power-law distributed. This 
method has been used to assess physiological neuronal avalanches and the results 
have shown that the power-law hypothesis is not rejected for this data [2]. It is there- 
fore of interest to determine whether this is also the case for the data from the model 
studied here. Briefly, this method finds the best fit to a power law of the distribution 
under study. The empirical data is then compared to distributions of the same size 
that are generated by randomly drawing values to follow the best-fit power-law dis- 
tribution. A p -value is calculated as the proportion of times that the empirical data is 
a better fit to the power law than the generated data (using the Kolmogorov-Smirnov 
test). As per Clauset et al. [1] the hypothesis (that the data comes from a power law) 
is rejected if the /?-value is less than 0.1. As we have observed (Figs. 4, 9), the dis- 
tribution of avalanche sizes appears to exhibit partial scale-free behaviour for low 
levels of external input (h = 0.1/ N, 0.01 /N) and the IAI distribution appears scale- 
free over a range of scales for h = 1/N. As in the companion paper [24], we fit a 
truncated power-law distribution up to an avalanche size of x max = j^N in the case 
of avalanche size distributions. We fit a power-law distribution without truncation to 
the IAI distribution. Testing the entire avalanche size distributions (consisting of over 
900,000 avalanches) yielded p = 0 indicating that the hypothesis that the distribu- 
tion follows a power law should be rejected. Similarly, taking the IAI distribution for 
h = 1/N, testing the whole distribution of over 6,000,000 IAIs (note that there are 
more avalanches and therefore IAIs with larger h due to the higher firing rate) yielded 
p = 0. Testing instead the first 100,000 avalanches (a similar order of magnitude to 
the number of neuronal avalanches tested experimentally) with h = 0.1/N yielded 
p = 0.46 indicating instead that the power-law hypothesis should not be rejected. 
Similarly, for h = 0.01/ N testing the first 10,000 avalanches yielded p = 0.13. These 
results are similar to those of the companion paper, where the power-law hypothesis 
was not rejected when the number of avalanches included in the distribution was of 
the same order as those tested experimentally, and they are indicative of the partial 
scale-free behaviour of the system in proximity to the critical regime. 

In the case of the IAI distribution testing the first 100,000 IAIs yielded p = 0.44 
indicating that a power law is a good fit to the data. Given that in this case we know 
that the IAI distribution is not a power law (and is in fact a weighted sum of hy- 
poexponentials), it is interesting to note that the hypothesis that the data follows a 
power law is not rejected when the number of data points is of the same order as that 
which have been tested experimentally, an observation that will be explained in the 
Discussion. When the power-law hypothesis is not rejected, Clauset et al. [1] employ 
a model selection process to determine the best model for the data. We did not carry 
out this testing here (as, at least in the case of the IAI distribution, we already know 
what the distribution is) and it may be that such a process would suggest that a power 
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law is not the best fit to the data. However, the results here (and those of the compan- 
ion paper) are indicative of the partial scale-free behaviour exhibited by the system 
in the region of the critical regime. 

2.5 Long-Range Temporal Correlations 

As discussed in the Introduction, long-range temporal correlations are another pos- 
sible signature of a system at (or near) a critical state and have also been observed 
in neurophysiological data [4, 1 1-15]. It is therefore of interest to determine whether 
this finite-size neuronal system with external input displays LRTCs — given that it is 
in the region of a critical regime — and whether LRTCs relate to other signatures of 
criticality, i.e. the presence of partial scale-free behaviour in the data distributions 
themselves. The latter is of particular interest given that we have seen a change in 
distributions as the system approaches the critical regime. As within any single simu- 
lation the level of external input is constant (and so does not itself display LRTCs) it is 
important to note at the outset that any LRTCs present in the dynamics of the system 
would be intrinsic to the system. Furthermore, the appearance of a power law within 
the distribution of any data set does not imply that the data will exhibit LRTCs and 
vice versa. (Consider points drawn at random from a power-law distribution — such a 
data set would not exhibit LRTCs.) 

In neurophysiological data, LRTCs have been observed in fluctuations of oscilla- 
tion amplitude (i.e. within continuous data) [4, 11-15] and also in discrete burst ac- 
tivity in our recent analysis of the inter-event intervals of bursts of nested oscillations 
in EEG recordings of extremely preterm human neonates [35]. Moreover, LRTCs in 
discrete data has previously been investigated by Peng et al. [44] and a number of 
other authors, for example [45-47], in their analysis of inter-heartbeat intervals. As 
the data from the model analysed here is discrete avalanche activity, we follow the 
approach of this previous analysis of LRTCs in discrete data, examining LRTCs in 
waiting times, i.e. in IAIs. 

We assessed the presence of LRTCs in IAIs through estimating the Hurst expo- 
nent, H , which describes the degree of self- similarity within the data. A Hurst ex- 
ponent of H = 0.5 indicates that there are no correlations in the data or short-range 
correlations only, for example a white noise process, whereas a Hurst exponent of 
0.5 < H < 1.0 indicates LRTCs in the data. Additionally, an exponent of 1 corre- 
sponds to l/f noise [44]. We estimated the Hurst exponent using detrended fluc- 
tuation analysis (DFA) — an approach that has been shown to produce more accu- 
rate estimates of the Hurst exponent than some other approaches [48] and has been 
used previously to assess the presence of LRTCs in neurophysiological data sets 
[4, 11, 15, 35]. DFA is a graphical method whereby the average root mean square 
fluctuations across a box size are compared across different box sizes and the gra- 
dient of the line of best-fit is the estimate of the Hurst exponent (for more detailed 
methodology see Peng et al. [44, 49]). We used a minimum box size of 5, with 50 box 
sizes linearly spaced on a logarithmic scale up to a maximum box size of 1/10 of the 
length of the IAI sequence [50]. Calculations were carried out using the MATLAB 
code of McSharry [51]. 

Figure 1 1 shows example DFA plots for IAIs from three simulations with a = 
w = 1, h = 1/N. It is important to notice from these plots that there is not a single 
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Fig. 11 DFA plot examining the presence of temporal correlations in IAIs. Plot of the average fluctuations 
F(n) against box size n for IAIs from simulations with a. = w = 1, h = 1/N and TV = 800 (a), TV = 3200 
(b) and TV = 172,800 (c). The data is best fit not by a single linear trend but by three lines (red, green, 
blue) between two crossover points (dashed black lines). For smaller box sizes the Hurst exponents (slope 
of the line — as annotated next to the individual lines: the DFA exponent a) indicates that correlations 
extend across these regions. However, for larger box sizes the exponents are closer to 0.5 suggesting 
that the correlations do not extend across these larger box sizes. With a larger system size (c) the upper 
crossover point increases to larger box sizes. With increasing system size the system approaches the critical 
regime—^ = -0.07 (a), X = -0.04 (b) and X = -0.005 (c) 



linear trend across all box sizes. Hu et al. [50] discussed the importance of identifying 
crossover points — box sizes at which there is a change in the linear fit of the data — 
within DFA plots. Failure to examine these trends leads to erroneous estimates of the 
Hurst exponents. A single linear fit across all the points would give an estimate of the 
Hurst exponent for that sequence. However, crossover points indicate that the same 
correlations (i.e. temporal behaviour) do not extend across the whole sequence. In the 
DFA plots here there are in fact three regions, each with a different linear trend, be- 
tween two crossover points. The best-fit to the data by three linear regions was found 
using the nonlinear regression function 'nlinfit' in the MATLAB environment, there- 
fore determining the crossover points. In Fig. 11(a), the Hurst exponent (slope of the 
line) of the first two regions (at smaller box sizes) are 0.83 and 0.62, respectively — 
exponents which indicate the presence of LRTCs within the data. However, the third 
region across the largest box sizes has an exponent of 0.51 which is close to the value 
of 0.5 which would indicate that there are no correlations in the data. This change in 
the exponents therefore suggests that the correlations shown in the data across small 
box sizes extend to a point but do not extend across the entire sequence length. 

When examining the presence of LRTCs it is standard practice to compare the 
exponent of the actual data to the exponent of the data randomly shuffled [4]. By 
shuffling the data this should destroy any correlations present and so the exponent of 
the shuffled data is expected to be approximately 0.5. We compared the original se- 
quence (whose DFA plot is shown in Fig. 11) with 500 shuffled sequences. The DFA 
plots for the shuffled sequences (data not shown) did not exhibit crossover points, 
with the same linear trend being observed across all box sizes. The mean exponent of 
the shuffled sequences was 0.50 with a range of 0.48-0.52. Therefore, as the expo- 
nents of the original sequence (at smaller box sizes) do not fall within the distribution 
of exponents for the shuffled sequences this further demonstrates that the original 
sequence exhibits complex temporal ordering with correlations that extend across a 
range of box sizes (up to the upper crossover). 
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Fig. 12 Changes with increasing system size, a The normalised (with respect to the largest box size) upper 
crossover box size increases with system size. The plot shows the average value across 10 simulations 
at the same system size and error bars indicate the standard deviation, b The IAI distributions and c the 
distributions of avalanche size for 6 different system sizes, scaled with respect to the mean IAI or avalanche 
size of each distribution, respectively. For all simulations a = w = 1, h = l/N 



2.5.1 Increasing the System Size 

As noted previously (see Fig. 2), with h = l/N as N -> oo the eigenvalue of the 
system A -> 0, i.e. the system approaches the critical regime with increasing system 
size. We might expect that as the system approaches the critical regime it is more 
likely to exhibit signatures of criticality and therefore that LRTCs would extend to 
larger box sizes as the system size is increased. We therefore investigated whether 
there was a change in the temporal correlations of the IAIs with system size, while 
maintaining all other parameters including h = l/N. For all system sizes investigated 
the DFA plots displayed three regions with different linear trends, as was discussed 
above. Figure 1 1 shows example DFA plots for the IAIs of three simulations with the 
smallest and largest system sizes examined. From this we can see that the pattern of 
the exponents in each of the cases remains the same — with the two lower regions hav- 
ing exponents indicative of LRTCs, while the exponent across the largest box sizes 
is closer to 0.5. Additionally, we find that with the larger system size the location of 
the upper crossover is at a larger box size. Figure 12 shows the level of the upper 
crossover (the crossover at the higher box regions) for different system sizes (from 
N = 3200 to N = 172,800) normalised with respect to the largest box size. The box 
size of the upper crossover increases with increasing system size. We did not find a 
change, other than small fluctuations, in the exponents for any of the three regions 
for different system size: across all system sizes, the mean exponent across the small- 
est box sizes (up to the first crossover point) was 0.73 with a range of 0.70-0.76. 
Between the first and second crossover points the average exponent across all simu- 
lations was 0.95 with a range of 0.89-0.99 (close to an exponent of 1 which would 
indicate l/f noise). The largest variation in exponents was for the region above the 
upper crossover point with an average exponent of 0.59 and a range of 0.46-0.73; 
on average this indicates that temporal correlations do not extend beyond the upper 
crossover. Thus, overall as the system size is increased the temporal correlations ex- 
tend across larger box sizes. This suggests that LRTCs will extend to infinite length 
(i.e. all possible box sizes) in the limit of system size — when the system reaches the 
critical regime. 
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Next we considered whether the distributions of IAIs and avalanche size them- 
selves changed with system size and whether the change in the correlation length 
observed above was reflected in a change in the distributions. Figure 12 also shows 
the IAI and avalanche size distributions for different network sizes. In both cases, 
the distributions for different system sizes have only small changes which can be 
accounted for by noise. Thus, as the system approaches the critical regime through 
increasing the system size there does not appear to be a change in the distributions de- 
spite the change in the temporal correlations. Moreover, LRTCs are present in the data 
but the distribution of avalanche sizes does not exhibit scale-free behaviour, i.e. these 
markers of criticality do not occur simultaneously in this case. By contrast, through 
approaching the critical regime by lowering the external input we have shown that the 
avalanches are more distinct and the distribution of avalanche sizes exhibits partial 
scale-free behaviour. 

2.5.2 The Effect on LRTCs of Decreasing the Level of the External Input 

We also examined the DFA exponents at the lower levels of external input (h = 0.l/N 
and h = 0.01/ N). In both cases there were no crossover points with a single lin- 
ear trend across all box sizes — data not shown. The exponents were 0.50 (range 
0.49-0.51, across 10 simulations with N = 800) and 0.56 (range 0.55-0.57) for 
h = 0.01/ N and h = 0.1/ N, respectively. Thus, at the lowest level of external input 
the IAIs do not exhibit LRTCs and there is a slight increase in the exponent as the ex- 
ternal input increases. This suggests that as the system approaches the critical regime 
through a decrease in the external input the temporal correlations are lost. Thus, the 
existence of LRTCs as the system approaches the critical regime is dependent on how 
the critical regime is approached, i.e. the region of parameter space — approaching 
the critical regime through increasing the system size extends the temporal correla- 
tions whereas decreasing the external input leads to a loss of long-range correlations. 
Moreover, this signature of criticality is independent from the other marker we have 
investigated — the presence of scale-free behaviour in the avalanche size distribution. 
Considering avalanche size and duration scale-free behaviour is present for the low- 
est level of external input, at which point LRTCs are lost. Thus, we find that markers 
of criticality are not only dependent on the region around the critical regime but also 
may not be present for the same parameter set. 



3 Discussion 

This paper specifically examined a driven finite- size neuronal system without an ex- 
plicit separation of timescales between the external input and the timescales of the 
avalanches themselves (i.e. while the external input we have considered was small, 
it was not limited to seeding the system at times of quiescence). By analytically tun- 
ing the system to be in the region of a critical regime we were able to examine the 
type of dynamics displayed by such a system and to investigate whether the dynamics 
displays signatures of criticality. In summary, we have shown that: 
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1 . As the system approaches the critical regime through a reduction in the external 
input the avalanches become more distinct and the distribution of avalanche sizes 
displays scale-free behaviour. 

2. With h = l/N the IAIs exhibit temporal correlations which extend across a range 
of bin sizes to an upper crossover. As the system approaches the critical regime 
through increasing the system size the length of the temporal correlations is ex- 
tended across a wider range of bin widths. These correlations (one noted signature 
of a critical system) are observed despite the fact that the distribution of avalanche 
sizes does not exhibit scale-free behaviour and does not change with the increase 
in system size. These temporal correlations are lost if the critical regime is instead 
approached through reducing the external input. 

3. The distribution of IAIs was theoretically derived and was shown to be a weighted 
sum of hypoexponentials. However, for h = 1/N (when the number of avalanches 
considered was of the same order as those tested experimentally) the hypothesis 
that the IAI distribution follows a power law was not rejected by statistical testing 
indicating the scale-free nature of the distribution at this level of the external input. 

3.1 Validity of the Model 

The model considered in this paper was a highly simplified neuronal system with a 
number of assumptions, such as equally weighted synapses and continuous constant 
external input. These assumptions were necessary in order to analytically tune the 
system to be in the region of a critical regime. Therefore, while this should not be 
taken as an accurate model of a neuronal system it is important that we first consider 
models such as this, examining markers of criticality, which will then aid our under- 
standing when building on this work with more complex models. This paper opens 
the way for future work examining the role of external input on signatures of critical- 
ity and the importance of the region of parameter space on network dynamics. Future 
work should also investigate the effect of topology on the dynamics [29, 52] and the 
effect of external input with different temporal and spatial characteristics. 

3.1.1 Purely Excitatory Synaptic Transmission 

The synaptic connections investigated in this model were purely excitatory. This not 
only simplifies the model for analytical investigations but is also of interest from a 
neurological perspective in terms of early brain development. Before birth, GABA 
is thought to have a depolarising effect on postsynaptic neurons and it is not until 
the nervous system reaches a more mature state that this neurotransmitter becomes 
inhibitory [53, 54]. While presynaptic inhibition is thought to be present at all de- 
velopmental stages [55] this effect can be considered to be taken into account in the 
model by the fact that neurons cannot re-fire until they have returned to the quies- 
cent state (i.e. inhibition in the model relates to the rate a at which neurons return to 
the quiescent state). We have recently shown that EEG recordings from very preterm 
infants (when GABA is still thought to be purely excitatory) exhibit LRTCs in the 
temporal occurrence of bursts of activity [35]. The model studied here may be a can- 
didate mechanism for the generation of this temporal patterning in the discontinuous 
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activity of the developing brain. Moreover, it is interesting to note that despite the 
fact that the system has purely excitatory postsynaptic connections and input, for 
these parameter regions, the model does not exhibit runaway excitation (saturation) 
but is able to maintain stable dynamics through the 'balance' of individual neuronal 
dynamics resulting from a trade-off between the rates at which neurons become ac- 
tive and quiescent. Indeed, while a number of authors have suggested that a balance 
of excitation and inhibition in neuronal networks leads to critical behaviour [56], the 
work here and in the companion paper shows that excitatory networks (i.e. networks 
without inhibitory neurons) can display the same behaviour. It can be speculated that 
this type of balanced activity in the region of a critical regime might be a way in 
which the brain avoids (for the most part) epileptic behaviour during early develop- 
ment, although it can also be argued that the decay rate contributes to a "balance 
condition" between excitation and inhibition [37]. 

3.1.2 The Activation Function 

Here we used a linear activation function for the transition of neurons from the qui- 
escent to the active states. However, physiologically neurons behave more like a sat- 
urating function. The linear activation function used here was chosen so as to be ana- 
lytically tractable and is also equivalent to a saturating function when input is small. 
However, considering instead a saturating function (see Appendix C) we found that 
the dynamics in the region of the critical regime shows similar behaviour to the sys- 
tem with the linear activation function. 

With both the linear and saturating activation functions, the critical regime can 
only be reached exactly in the absence of external input. A positive external input 
therefore drives the system away from the critical regime. However, with a quadratic 
activation function (see Appendix C) the system, with a positive external input, has a 
critical fixed point and the system can be tuned directly to this regime. With this acti- 
vation function the dynamics does not appear to exhibit burst like behaviour, however, 
analysis shows that the activity fluctuates about the critical regime in an 'avalanche- 
like' manner. Thus, while a quadratic function does not best describe activation in 
a neuronal network, we may further conclude that signatures of criticality are not 
universal and can be examined only in relation to the specific critical regime of the 
system (see Appendix C). 

3.1.3 The Binning Approach 

As described previously, the binning method separated avalanches where the time 
difference between consecutive spikes was greater than the average time difference 
between consecutive spikes across the entire simulation. This was the approach taken 
by Benayoun et al. [38]. However, it is worth noting that this is a slightly differ- 
ent approach to the method that has been used experimentally to separate neuronal 
avalanches — first proposed by Beggs and Plenz [5, 6]. In their analysis neuronal firing 
is distributed into bins of width of the average time difference between consecutive 
spikes (8t) and firing is separated into avalanches by bins in which no firing occurs. 
Thus, two spikes may be greater than the average time difference 8t apart but still 
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remain in the same avalanche if they fall within consecutive bins. The theoretical 
derivation of the IAI distribution relied on the fact that all consecutive active to qui- 
escent transitions or single quiescent to active transitions with a length greater than 
the average time between two spikes is an IAI. This would not be the case if the 
alternative (Beggs and Plenz) binning approach was used to determine avalanches. 
If this alternative approach had been used the distributions of consecutive active to 
quiescent transitions and single quiescent to active transitions would be the same, 
but transitions of length slightly greater than or equal to the average time between 
consecutive spikes (in fact up to twice this average) may or may not form part of the 
IAI distribution depending on the exact binning. It is also important to note that with 
the binning method used here, even with dense neuronal firing (which occurs if the 
external input is increased from the levels studied here), as there is always an average 
time between consecutive spikes it is always possible to separate the dynamics into 
'avalanches'. 

Additionally, both these binning approaches differ from that used in non-driven 
systems such as the classical sand-pile model [18] and the system investigated in the 
companion paper to this article [24]. In those models an avalanche consists of all 
firing until the system returns to the fully quiescent state and so, for example, the 
system may have a long period without firing in which neurons switch to the inactive 
state but this will not be designated as two separate avalanches (if the system has 
not returned to the fully quiescent state) even when the period exceeds the average 
difference between consecutive spikes. Future work is needed to fully investigate how 
the differences in these avalanche definitions affect the distributions of size, duration 
and IAIs and care needs to be taken when interpreting the results from these different 
approaches. 

3.1.4 Validity ofDFA and the Investigation ofLRTCs 

DFA is one method by which to estimate the Hurst exponent and was chosen here as it 
has been shown to be an accurate estimate [48]. Moreover, it is a graphical approach 
and so can be used to check for crossover points [50]. As the Hurst exponent can 
only be estimated it is considered to be best practice to check the consistency of the 
exponents using two methods [57]. However, as non-graphical methods only give 
single numerical values they cannot be interpreted when crossover behaviour exists. 
Given that there were crossover points we only considered DFA with this analysis. 

Crossover points within a DFA plot have been shown to exist when the same corre- 
lations do not extend across the whole data sequence in analytically constructed data 
[50]. The crossover points in the data here can therefore be interpreted in this way, as 
points at which the correlations in the sequences change. It is important to understand 
that these crossover points (and box sizes in general) relate to the sequence length. 
For example, a box size of 10 indicates detrending across 10 consecutive IAIs. As the 
IAIs themselves can be of variable length the box size does not relate to a particular 
simulation time. Future investigation is needed to determine the relationship between 
the model and crossover points. 

Correlations extended across a range of box sizes with this range extending as the 
system size increased and the system approached the critical regime. It appears that 
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correlations would extend across an infinite box size in the limit of system size. Thus, 
as the critical regime is approached in this way, this signature of a critical system 
emerges. LRTCs have been demonstrated previously in discrete neurophysiological 
data, in the waiting times of burst activity in cultures [34] and in the bursts of ac- 
tivity recorded using EEG in very preterm human neonates [35]. To our knowledge, 
waiting times of neuronal avalanches have yet to be examined in this way. How- 
ever, such a study would provide an additional link between studies on the neuronal 
scale and studies on a wider network scale for which LRTCs have been observed 
in the fluctuations of oscillation amplitude. Palva et al. demonstrated strong correla- 
tions between power-law exponents of avalanche size distributions and exponents of 
LRTCs in fluctuations of oscillation amplitude in human MEG recordings [33]. Re- 
cent computational work also demonstrated a link between neuronal avalanches on 
the one scale and LRTCs on a wider temporal scale and the authors called for future 
work in this area [32]. However, the authors of this study did not investigate LRTCs 
in the waiting times of the avalanches themselves. Interestingly, in the model studied 
here, LRTCs were observed when h = 1/N but not for lower levels of external input. 
Thus, they were not observed when the avalanche size distribution exhibited scale- 
free behaviour — the type of distribution observed for avalanches recorded in vivo and 
in vitro [5, 6, 9]. It would therefore also be interesting to assess whether altering the 
driving force experimentally in vitro would lead to the types of dynamics (LRTCs) 
observed here. 

3.2 Partial Scale-Free Behaviour in Avalanche Size 

Statistical testing of the avalanche size distribution (with h = 0.l/N, 0.01 / N) did not 
reject the hypothesis that the distribution followed a power law when the number of 
points within the distribution was of the order of the number of avalanches recorded in 
the experimental setting. Only with larger numbers of avalanches was the hypothesis 
that the distribution is a power-law rejected. This is to be expected — as has been 
discussed by Klaus and Plenz [2], when a distribution deviates from the expected 
distribution by more than noise from sampling then given a large enough number 
of samples the power-law hypothesis will eventually be rejected. The fact that the 
power-law hypothesis was not rejected for lower numbers of avalanches demonstrates 
the partial scale-free behaviour of the system in the region of the critical regime. 
Moreover, this highlights the fact that stringent statistical testing, such as this, with 
high sampling may lead to rejecting the power-law hypothesis and so rejecting the 
criticality hypothesis even when the system is critical. 

3.3 Waiting Times 

In addition to increasing the physiological realism of the model, investigating the 
driven system also has the advantage of producing waiting times (in this case 
termed IAIs). In the companion paper the simple reseeding of the network with a neu- 
ron set to the active state implied that there was no waiting times between avalanches. 
Other authors have reseeded by increasing the membrane potential but stipulated that 
neurons must reach a threshold for them to become active (and a new avalanche to 
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start) [18, 25]. This does lead to waiting times, however, these are not the same as the 
waiting times investigated in this model which are intrinsic to the network dynamics 
rather than as a result of network reseeding. 

Recent work by Lombardi et al. [31] showed that the waiting times between 
neuronal avalanches recorded in vitro have a distribution with an initial power-law 
regime. The authors suggest that the shape of the distribution relates to up and down 
states within the network (which exhibit critical and subcritical dynamics, respec- 
tively) and are able to reproduce the non-monotonic waiting time distribution in a 
computational model in which neurons switch between up and down states depend- 
ing on short-term firing history. Interestingly, the distribution they observe is similar 
to the IAI distribution for the system with h = 0.1/ N, see Fig. 9(b), which also has a 
scale-free initial regime albeit over a shorter range to that presented by Lombardi et 
al. It is therefore possible that the waiting time distribution observed experimentally 
fits with the model constructed here. It would be interesting to investigate whether a 
change in input to the network in vitro alters the distribution in a similar way to those 
distributions seen in Fig. 9. 

Additionally, for different parameter ranges different distributions were observed, 
in the IAI distribution as well as the distributions of avalanche size and duration. This 
leads us to the important conclusion that power-law distributions will not necessarily 
be displayed by systems in the region of a critical regime. Therefore, this work sug- 
gests that the absence of a power law in experimental data should not necessarily be 
taken to conclude that the system does not lie in the region of a critical regime. This 
was also seen in the companion paper where it was shown that despite being analyti- 
cally tuned to the critical state (without the presence of external input) the avalanche 
size distribution was not a power law although it did exhibit partial scale-free be- 
haviour. The fact that the system may not exhibit power laws when close to (or at) the 
critical regime is an important finding given that the system is of finite size as will be 
the case in the experimental setting. This highlights the necessity of examining other 
markers of criticality before conclusions about the critical nature of a system can be 
drawn. 

3.4 Dynamic Range and Power Laws 

Coinciding with results from previous authors [16, 29] we showed that the system 
exhibits optimal dynamic range when the branching parameter is equal to one. When 
calculating the dynamic range of a system, we emphasised that this value was de- 
pendent on the critical state of the system calculated when there was zero external 
input. We have shown that tuning a system to this critical point but then driving it 
with different levels of external input has considerable effect on the distribution of 
avalanche sizes. For non-zero h the corresponding ODE would, in the strictest sense, 
not be considered critical. Importantly, however, tuning to the critical point of the 
system with zero external input, maximises the dynamic range. 

Dehghani et al. [58] showed that in vivo (contradictory to the results of Petermann 
et al. [9] and Hahn et al. [8]) avalanches were not well approximated by power laws, 
but they were more likely to approach exponential distributions. They contrast this 
with the evidence that the brain is operating at criticality from in vitro studies [5, 59] 
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where avalanches are well approximated by power laws. Here we argue that external 
input and functional benefits [17] such as dynamic range, information transmission 
and information capacity, provide an interesting possibility as to the reason why in 
vivo and in vitro studies could potentially give different results. The critical brain hy- 
pothesis demands that in isolation from its natural surroundings (in vitro) and whilst 
having no external influences acting upon it (akin to the model with h = 0 we studied 
in the companion paper [24]), a culture should exhibit signs that it is tuned to crit- 
icality (i.e. avalanches that are well approximated by power laws). However, when 
observed in vivo, and thus with external inputs acting upon it, a critical brain may no 
longer exhibit avalanches approximated by power laws but would instead optimise 
functional benefits such as the dynamic range and information transmission [17]. In 
our model we have shown that tuning the parameters to the critical regime does in- 
deed maximise the dynamic range, but it is the level of external input that dictates 
whether the avalanche distributions exhibit partial scale-free behaviour. For this rea- 
son, avalanches recorded in vivo that lacked a power-law distribution would not be 
contradictory to criticality but instead an expected result. This further supports our 
suggestion in the companion paper [24] that future work should shift focus away 
from characterising avalanche distributions to more appropriate metrics. 

3.5 Two Routes to Criticality 

In this paper we examined two different parameter changes such that the system ap- 
proaches the critical state: increasing the system size and lowering the overall level of 
the external input. Despite the fact that in both cases the critical regime is approached, 
the dynamics and the signatures of criticality observed are different. With increasing 
system size the temporal correlations extend across a wider range. However, the dis- 
tributions of the avalanche characteristics remain the same and the distribution of 
avalanche size does not exhibit scale-free behaviour. By contrast, for lower overall 
levels of the external input the distributions of avalanche size and duration do exhibit 
partial scale-free behaviour. However, in this case as the critical regime is approached 
the temporal correlations in the avalanches are lost. At these lower levels of the exter- 
nal input we also observe a greater separation of the avalanches suggesting that the 
avalanches have less of an influence on each other which would explain this loss of 
LRTCs. Thus, as the system approaches the critical state in two different regions of 
the parameter space the dynamical properties of the system are very different. Sig- 
nificantly, this implies that not just the critical state alone but the region around the 
critical regime is an important factor in the system's dynamics. 

In conclusion, we have shown here and in the companion paper that in a finite- 
size neuronal system in the region of a critical regime the distributions of avalanche 
attributes need not be a power law. The current assumption in the literature is that 
power-law dynamics implies criticality and vice versa that systems without power- 
law dynamics are not in the region of a critical regime, however, the results here 
suggest that this assumption need not be true. Moreover, we found that long-range 
temporal correlations and scale-free distributions are not dependent on proximity to 
the critical regime alone but on the region of the parameter space. The results further 
highlight the need for future work examining the type of dynamics we might expect 
from such systems. 
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Appendix A: Dynamic Range 

Whilst [16] and [29] consider a discrete model where multiple events can happen per 
time step, here we show analytically that our continuous model will exhibit the same 
maximisation of the dynamic range when Ro = 1 . Here we use the calculation of Ro 
for a system where there is no external input (h = 0) and thus Rq = w/ct. 

We begin by defining (as in Kinouchi and Copelli [16]) F mdx (Ro) as the satura- 
tion level of neurons in a network assuming a large external input h. For our model 
^max(^o) = N for all Ro. Similarly we define Fo(Ro) as the steady state solution of 
the mean field ODE for the system when there is zero external input, i.e. 



Additionally let the response function (approximating the mean firing rate) F x (Ro) = 
Fo(Ro) +x[F max (Ro) ~ Fo(Ro)] [16] giving 




(N — A) — aA. 

N 



Therefore, solving this we have 



Wo) = 



0 if#o<l, 
if*o>l. 



F x (Ro) = 



Nx 



if Ro < 1, 



N[l 



a 



(1-jc)] if/? 0 >l. 



w 
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Finally, let A(a, y) be the number of active neurons at the steady state in a regime 
where Ro = o and h = y (where o and y are dummy variables and h is the external 
input), then the dynamic range A(Rq) is defined (similarly to [16]) as 



A(R 0 ) = 




ho.i is the level of external input such that A(Ro, /zo.i) = Fo.i(^o) = ^b.l and 
/*0.9 is the level of the external input such that A(Ro, ho.9) = Fo.9(Ro) = F0.9. 

We note that in [16, 29], the logarithm of this is taken but as the logarithm is an 
increasing function it is unnecessary to scale in this way for the result we obtain. 
Whilst using Fo.i and F0.9 is the standard for calculating the dynamic range these 
values are somewhat arbitrary [16] and can be generalised to k\ and k 2 , respectively. 
To calculate the dynamic range analytically we consider the two regimes of Ro, firstly 
Ro < 1 and secondly Ro > 1 . 

*o<l 

Here the steady state is given by 



N 



+ h k j(N-F k )-aF k = 0 

(wk + h k )(N - Nk) - aNk = 0 
ak 



hk 



l-k 



— wk, 



thus 



hk 2 
hki 



[ak2 - wk2(l - &2)](1 - k\) 



[ak\ — wk\(l — k\)](l — £2) 

k 2 (l-ki)[l-Ro(l-k 2 )] 
k 1 (l-k 2 )[l-Ro(l-h)Y 



Here the steady state is given by 



ujFk_ 
N 



+ h k )(N-F k )-aF k = 0 



w[ 1 - -(1 -k) 







~Na 






) +h k 




— (l-k) 


-aN 
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hk 
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(w — a + ak) 
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thus 

_ hk 2 _ &2(1 - a + afo) 

hk x k\(\ — k2)(w — a + ak\) 

= k 2 (l-k 1 )(Ro-l+k 2 ) 
~ ki(l-k 2 )(Ro-l+h) 

A.l Maximum of A (Rq) 

Calculating the derivative of A(Rq) we find that if 0 < k\ < k 2 < 1, then for Ro < 1, 
> 0, whilst for 7?o > 1, ^ < 0- Thus, there is a critical point at Ro = 1 where 
the maximum of A(Ro) is achieved — see Fig. 1. It is worth noting that A(Rq) is 
independent of N and only depends on the choice of k\ and k 2 . 



Appendix B: Driving the System from a Subcritical and Supercritical State 

Throughout the paper we have examined parameters such that the system is critical 
when there is no external input. In the presence of a small external input we therefore 
investigate driving the system in the region of this critical state. In the companion 
paper [24], with no external input, we also investigated the system with subcritical 
and supercritical parameters. In this appendix we briefly examine the dynamics of the 
system as it is driven from these states by an external input. 

Figure 13 shows raster plots of network firings when the system is driven from 
a subcritical and supercritical state with h = \/N . Compared with the critical case, 
see Fig. 3(a), with the subcritical parameter set the bursts appear to be shorter and 
consist of fewer neurons firing. Conversely, in the supercritical case the bursts appear 
longer and consist of denser network firing. Figure 14 shows the IAI distributions 
for the subcritical and supercritical parameters. As expected from the raster plots, the 
IAIs are longer in the subcritical case compared with the critical (Fig. 9(a)) and the 
supercritical. While the subcritical distribution appears to exhibit partial scale-free 
behaviour similar to the critical case, the supercritical distribution loses this appear- 
ance. The distributions from simulations are shown with the theoretical distribution 
calculated as previously described as a weighted sum of hypoexponentials. 

Figure 15 shows the distributions of avalanche size and duration in the subcritical 
and supercritical cases. Contradicting what we would expect from the raster plots 
we find that the avalanche sizes are smaller (on average) in the supercritical system. 
In the companion paper we showed that the supercritical distribution (without the 
presence of external input) had an increased number of large avalanches compared 
with the distribution for the system at criticality. However, we do not find this here. 
As the firing with the supercritical parameters is relatively dense we believe that this 
highlights a limitation with the binning method in this case. We suggest that future 
research should focus on how binning can influence avalanche distributions. 
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Fig. 13 Raster plots of neuronal firing for the network driven from subcritical and supercritical states. The 
network firing for subcritical, a = 1.1, w = 1 X < 0 (a) and supercritical, a = 0.9, w = 1 X > 0 (b) 
parameter sets. Here we investigate the system with a small external input (h = 1/N) which drives the 
system slightly away from these fixed points. The red line indicates the level of firing in 1 ms bins. The 
subcritical case appears to give rise to smaller bursts and the supercritical case leads to a greater level of 
firing and longer burst activity compared with the critical system (see Fig. 3) 




Fig. 14 Distribution of IAIs for the system driven from subcritical and supercritical states. The theoretical 
(red) and simulated (black) IAI distributions with a = 1.1 (subcritical) (a) and a = 0.9 (supercritical) (b) 
parameters. The blue line in a indicates a linear fit, i.e. a fitted power law with an exponent of 2.37. These 
distributions are with TV = 800, w = 1, h = 1/N and the theoretical distributions were calculated up to a 
level of initial active neurons which occur with a cumulative probability of 0.9 and 0.13, respectively (see 
main text) 



Fig. 15 Distributions of 
avalanche size and duration for 
the system driven from 
subcritical and supercritical 
regimes. Avalanche size and 
duration distributions for the 
system driven from subcritical 
(a = 1.1) (a, c) and supercritical 
(a = 0.9) (b, d) states. 
Simulations were run with 
TV = 800, w = 1, h = 1/N. The 
red line in a shows the linear fit 
with a slope of -2.66 
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Fig. 16 Raster plots for different levels of the external input with the saturating activation function. Neu- 
ronal firing during the first 2 seconds of example simulations with h = 1/N (a), h = 0.1 /N (b) and 
h = 0.01/ N (c). For all three simulations w = 1, a = 0.25 and TV = 800. Comparing with Fig. 3 we find 
that while in this case the firing rate is lower (note the longer time scale over which the raster plot is dis- 
played) the overall pattern is the same, with the avalanches becoming more distinct with lower levels of 
the external input 



Appendix C: Altering the Activation Function 

Throughout this paper we considered a linear activation function. What happens if a 
different activation function is chosen? Do we observe the same type of dynamics? 
In this appendix we briefly investigate two other activation functions: an exponential 
and a quadratic. 

First let us consider the system with an exponential activation function such that 



dA ( 1 1 



(N -A)-aA. 



This function saturates and so is somewhat more realistic than the linear function 
considered previously. Also, note that the function is set such that when A and h 
are both zero we also have f(x) = 0, i.e. without any external input and with no 
active neurons the network will remain in this state. With this activation function the 
eigenvalues of the fixed points are given by 



X = f(x)(N-A)-f(x) 



We have 



{w/N)e- ((w > N)A+h) w( 1 \ / 1 \ w/l , 

f <*> = (1 + e -((uVAW))2 = n { fix) + 2 J \2 ~ fiX) ) = N \4 ~ f {X) 

and so this could be used to find a critical fixed point along with the fact that at the 
fixed point of the system we have 

aA 

/(*) = 



(N-AY 



which defines the level of the external input at the critical fixed point. 

As before, consider initially the case where there is no external input (h = 0). In 
this case A = 0 is a fixed point, which is critical (with X = 0) if and only if a = w/4 
by the above equations. What happens to this system in the presence of small exter- 
nal input? Figure 16 shows the raster plot for the three different levels of the external 
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Fig. 17 IAI, avalanche size and duration distributions for the system with the saturating activation func- 
tion. The distributions are for simulations with h = l/N (a, d, g), h = 0.1/ N (b, e, h) and h = 0.01/ TV 
(c, f, i). For all simulations TV = 800, a = 0.25, w = 1 and the distributions are pooled from 10 simu- 
lations each of length 10 4 seconds. The red lines indicate linear fits on the double logarithmic scale, i.e. 
fitted power laws with exponents of 2.65 (a), 1.81 (e), 1.49 (f) and 1.84 (h) 



input considered previously: h = 0.01/ N, 0.1/ N, l/N. Comparing with Fig. 3, the 
firing rate is lower with the saturating function studied here, however, the overall 
pattern of firing is the same. For all three levels of the external input we continue 
to observe avalanche dynamics and for lower levels of the external input (as the sys- 
tem approaches the critical regime) these avalanches become more distinct. Figure 17 
shows the avalanche size, duration and IAI distributions for each of these three levels 
of the external input. Comparing with Figs. 4 and 9 we find that a similar relationship 
with the critical regime emerges. With h = l/N the IAI distribution shows scale-free 
behaviour (note that by the same derivation as previously, theoretically the distribu- 
tion is a weighted sum of hypoexponentials). For lower levels of the external input 
the scale-free behaviour in the IAI distribution is lost but the distribution of avalanche 
sizes appears scale-free. As was shown previously for the system with a linear acti- 
vation function, we also found that when h = l/N the IAIs exhibited LRTCs up to a 
crossover point (data not shown). For lower levels of the external input these correla- 
tions were lost. 

With both the linear and saturating activation function we considered the system 
in the region of the critical regime, with the system driven from the critical regime 
by the positive external input. Consider the system instead with a quadratic activation 
function: 

dA ( w 9 \ 

~d7 = \N + ) (N ~ A) ~ aA ' 

Springer 



Page 38 of 42 



C. Hartley et al. 



With this activation function the fixed points are given by 




(h+a)A + hN = 0, 



with eigenvalues 



A = g\A) = -3 — A 2 + 2wA -h-a. 



Solving these simultaneously we find 




which define the parameter space and the value of the fixed point for which a critical 
fixed point can be obtained. Thus, we find that unlike the model with the linear (and 
saturating) activation function, here with a non-zero external input it is possible to 
tune the system so that it is directly at the critical regime. 

Upon examining this parameter space one can note that in many cases there also 
exists a stable (positive) fixed point as well as the critical fixed point. From simulating 
such a system we found (data not shown) that the dynamics of the system is quickly 
attracted to the stable fixed point and so the critical fixed point has little affect on the 
dynamics. Therefore, to have a system which is affected by a critical fixed point in the 
presence of a non-zero external input (in the case of this activation function and where 
positive parameters are required) the critical regime must be the only fixed point of 
the system. Given that g(A) is a cubic equation, to achieve a single fixed point which 
is critical this point must be an inflection point with g'(A) = 0 and g"(A) = 0. From 
these equalities we find that the critical fixed point is A = N/3 and we must also have 
/* = ^anda = ^f. 

Figure 1 8 shows a raster firing plot and the number of active neurons throughout 
a simulation for the system with a single critical fixed point. As would be expected, 
the number of active neurons fluctuates about the critical point. Previously when 
considering avalanche dynamics we have binned the firing. However, as noted in 
the Discussion the binning method will always separate firing into avalanches and as 
there are no clear periods of inactivity this does not seem appropriate here. Recall that 
in the zero input case (see the companion paper) we seeded the system so perturbing it 
away from the fully quiescent state (which was the critical fixed point) and defined an 
avalanche as the firing that occurred before the system returned to the fully quiescent 
state. In a similar approach here it is possible to define an avalanche as the number 
of neurons that fire in a single excursion from the critical fixed point. We therefore 
counted the number of neurons that fired from when the system was deflected (either 
in a positive or negative direction) from the fixed point (A = N/3) until the next time 
at which the system had exactly N /3 active neurons. Figure 18 shows the probability 
distribution of the size of the avalanches defined in this way. The distribution appears 
to be scale-free over a range of scales. 

Thus, while critical dynamics may not be apparent initially when examining data 
(for example if we were to look at the overall dynamics from the simulations with 
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Fig. 18 Dynamics of the network with a quadratic activation function, a Raster plot of network firing 
for simulation with TV = 800, h = a = w = 0.01. b For the same simulation this plot shows 

the number of active neurons (A) at each simulation step. The distribution of avalanches (c) and positive 
avalanches only (d) — as described in the main text — pooled from 20 simulations of 1000 seconds in length. 
The red lines indicate linear fits, i.e. fitted power laws, both with exponents of 1.48 

quadratic activation function), we can observe signatures of criticality when the dy- 
namics is examined in relation to the known critical regime. Here we can note that 
the network firing fluctuates about the critical regime — that is, the number of active 
neurons fluctuates about this regime and so the average number of active neurons 
across the course of a simulation is approximately equal to the critical state of N/3. 
It might therefore be interesting to examine the fluctuations about the mean activity 
level in experimental settings where activity is continuous (i.e. cannot be described 
as intermittent avalanche-like activity) to determine whether signatures of criticality 
are present. Indeed, such an approach has been taken previously to examine MEG 
data, thresholding at the median level [60] . 
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